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ABSTRACT 

We discuss the use of recursive enumeration schemes to obtain low and high tem- 
perature series expansions for discrete statistical systems. Using linear combinations of 
generalized helical lattices, the method is competitive with diagramatic approaches and is 
easily generalizable. We illustrate the approach using the Ising model and generate low 
temperature series in up to five dimensions and high temperature series in three dimen- 
sions. The method is general and can be applied to any discrete model. We describe how 
it would work for Potts models. 
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Introduction 

Expansions about either infinite or vanishing coupling have long been a major tech- 
nique for the study of statistical systems and field theories. These series usually involve a 
diagrammatic analysis which becomes rapidly more complex as the order increases. Thus it 
would be interesting to have a fully automated technique for the generation of the relevant 
terms. 

Here we discuss a purely mechanical method to generate the low and high temperature 
expansions for discrete systems. The approach does not involve explicit graphs, but rather 
relies on a recursive computer enumeration of configurations. We illustrate the approach 
on Potts and Ising models, although it is considerably more general. 

The method is based on a recursive transfer matrix procedure of Binder [1] for the 
explicit solution of discrete models on small lattices. Enting [2] discussed how to combine 
such solutions on small lattices to obtain low temperature series. Guttmann and Enting 
have pushed this finite lattice method to obtain rather high order low temperature series 
for the three dimensional Ising model [3]. Our approach is similar in spirit to this work, 
although it differs in many technical details. In Ref. [4] these ideas were further developed 
in the context of finite size scaling and the analytic structure of the partition function. 
Ref. [5] explored using these exact counting procedures on helical lattices to extract the 
low temperature series. Helical lattices have been further generalized in Ref. [6] enabling 
one to calculate the low temperature series for the three dimensional Ising model to 50 
excited bonds. Ref. [7] applied these methods to Potts models in two and three dimensions. 

This paper is primarily intended to explain these methods in more detail and explore 
extensions. In a recent preprint Vohwinkle [8] has adapted the diagrammatic shadow 
method to obtain Ising and Potts expansions to several more terms than we have been 
able to obtain. Our method is, however, quite general and easy to implement. It is an 
open question whether some of the ideas of Ref. [8] can be adapted into our scheme to get 
even longer series. 

Recursive counting 

We begin with a discussion of the recursive approach to solving small systems exactly. 
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This section will also serve to establish notational conventions. We illustrate the basic 
method with the Ising model on a finite three dimensional simple cubic lattice. On each 
site % is a spin &{ taking the values ±1. The energy of the system is 

E=Y J Q-Wi) (!) 

where the sum is over all nearest neighbor pairs of spins, each pair being counted once. At 
inverse temperature (3, the partition function is the sum of the Boltzmann weight over all 
configurations 

Z = Y j e~" E (2) 
M 

Organizing the set of configurations by their energy, we rewrite this as a sum over E. This 
introduces the density of states function P(E) representing the number of states of the 
system with the given energy E. Thus, we have 

67V 

Z = P(E)u E/2 (3) 

where N is the number of sites and u = e~ 2 @ . If we consider, for example, an iV 3 lattice, 
there will be 2 nS states, but the solution for the partition function can be expressed in 
terms of 0(N 3 ) integers P(E). 

For a given lattice we compute the coefficients P(E) exactly using a transfer matrix 
to assemble the system one site at a time. This recursive construction enables us to 
build up a lattice with arbitrary length in one of the three dimensions. For the series 
analysis it is important to always continue the recursion sufficiently to avoid finite size 
errors in this "longitudinal" direction. At intermediate times the process requires an 
explicit enumeration of any exposed two dimensional slice. This effectively reduces the 
computational complexity to that of a system of one less dimension. Thus the solution of 
an A 3 lattice requires, at most, the explicit enumeration of only 2^ states. This enables 
us to work with sizes which would be impractical for an explicit enumeration of all states. 

The starting point of the method is a list of all states and corresponding energies 
for a single transverse layer of the lattice. All spins outside this layer are frozen to the 
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same value; that is, the boundary conditions in the longitudinal direction are cold. Spins 
are then sequentially freed to build up the lattice in this third direction. At intermediate 
stages the computer stores the exact number of states of any given energy and specified 
exposed top layer. Storing the top layer as the bits of an integer /, define p(E, I) to be 
this count. When a new spin or set of spins is added, we obtain the new counts p'(E, I) 
as a sum over the old counts 

p'(e,i) = J2p(e-A(i, I'), I'). (4) 
I' 

Here /' runs over all integers differing from / at most in the bits representing the newly 
covered spins, and A(I, I') is the change in energy from any newly changed bonds. In 
Ref. [5] the spins were added one layer at a time, while here we add them one at a time. 
Thus for the present calculation the sum in the above equation is only over two terms, 
representing the two possible values for the newly covered spin. After the lattice is grown, 
a sum over the top layers gives the resulting P(E) for the entire system 

P{E) = Y J P{EJ). (5) 
i 

The low temperature series 

Note that as the temperature goes to zero so does the variable u. Thus Eq. (3) is 
itself the low temperature expansion for Z . From it, we compute the corresponding series 
for the average energy per site, 

{E) = E^W = 2 » ) log(z) (6) 

Comparing this expectation value before and after the last spin is added, we obtain the 
increase in the average energy per new site. Expanding this in powers of u gives 

(E/N) = J2e jU i (7) 

3 

We are interested in the coefficients ej in the infinite volume limit. 

At zero temperature (f3 = oo) the only states which survive have all spins parallel. As 
the temperature increases, groups of spins can flip in this uniform background. A single 
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flipped spin has 6 excited bonds, and thus represents the first non-trivial term in our 
expansion. In our units, each excited bond has energy 2 and there are 6 such bonds for a 
single flipped spin; thus we have e§ = 12. Continuing to more complicated combinations 
of flipped spins gives the usual diagrammatic method to obtain the further coefficients. 

Note that any enclosed group of flipped spins must always have an even number 
of excited bonds. Thus the expansion only involves even powers of u. Our method of 
construction is such that when a spin is added to the helix, we account for the energy of 
both the forward and backward bonds at once. This, combined with the cold boundary 
conditions at both ends, ensures that we generate only even powers of u in our expansions. 

Helical lattices 

Computing Z exactly on a periodic lattice of size N x N x N, the order to which the 
weak coupling expansion for (E/N) will agree with the infinite volume limit is AN — 2. 
At this order a line of N flipped spins can wrap around the lattice and show finite size 
effects. Such a configuration will have energy AN rather than the AN + 2 it would have in 
infinite space. This is the smallest excitation affected by boundary effects and hence, on a 
periodic lattice of size N, the expansion is valid through 0(AN — 2). 

The order to which the series is correct can be increased by changing the boundary 
conditions to require more spins to be flipped to wrap around the lattice. Ref. [5] showed 
how a version of helical boundaries allowed an N x N transverse slice to be mimicked with 
only [(N 2 + l)/2] sites. Here we extend this idea to include the helicity into the direction 
in which the lattice is grown. 

We build our lattices one site at a time; so, it is natural to imagine the sites lying 
in a line. We do not, however, consider sequential sites as nearest neighbors. Instead, 
we introduce three integer parameters {h x , h y , h z } representing the distance along the line 
to the nearest neighbor in the corresponding x, y, or z direction. Labeling sites in the 
sequence by their ordinal number i, the nearest neighbors of site i are at i ± h x , % ± h y and 
i ± h z . For convenience, we assume 

h x <h y < h z (8) 
With this convention, as we grow our lattice, all sites more than h z steps back in the chain 
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are covered. Thus the recursive methods of Refs. (1-3) only require us to keep explicit 
track of the h z "exposed" spins at the end of our chain. The computational work grows 
exponentially with this number; thus, we wish to keep h z as small as possible. 

A minimal closed loop on such a lattice consists of a number of steps in each spatial 
direction such that 

n x h x + riyhy + n z h z = 0. (9) 

Here rii represents the number of steps in the i— th direction. The length n of such a loop 
is 

n = \n x \ + \n y \ + \n z \ (10) 

n is the "effective" periodic size of the lattice for our series construction, and, as argued 
above, the series is correct to 0(u 4n ~ 2 ). On an infinite cubic lattice the only solution 
to Eq. (9) is the trivial case n; = 0. On a finite lattice, any other solution represents 
a finite size correction. Flipping a chain of spins along such a closed path generates a 
state with 4n excited bonds, and creates a potential error in the series at that order. As a 
simple example, (h x ,h y ,h z ) = (3,4,5) with (n x ,n y ,n 2 ) = (1,-2,1) gives a minimal loop 
of length 4. Such a lattice will give the series equivalent to that on a 4 3 lattice, but with 
only h z = 5 sites in cross section. Similarly, (h x ,h y ,h z ) = (19,21,24) has closed loops of 
length 10 corresponding to (n x , n y , n z ) = (3, —5, 2). Here 24 sites mimic a 10 by 10 cross 
sectional lattice, thus saving a factor of 2 76 in computational effort. 

Note that Eq. (9) tells us that if we regard n and h as vectors, they are orthogonal. 
Thus a simple way to visualize our lattice is as an infinite one with all sites which lie in 
any single plane orthogonal to h as identified with each other. Fig. (1) attempts to show 
this construction. Considering the plane through the origin, all the sites lying in this plane 
themselves form a lattice. Closed loops that contribute finite size corrections consist of 
sets of flipped spins connecting the sites of this lattice. 

Cancelling loops 

We now discuss how forming linear combinations of the energy series coefficients from 
a set of finite helical lattices can give the infinite volume series to a higher order than any 
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individual lattice in the set. The approach here differs in details but is similar in essence 
to the combining of partition functions in the finite lattice method of Refs. [2] and [3] . 

Given a set of parameters (h x , h y , h z ), it is straightforward to enumerate the minimal 
closed paths. A different set of parameters corresponds to a different set of such paths. 
However, any erroneous contribution to the coefficients from a particular such path is, 
by symmetry, independent of any permutation or sign changes in the numbers (n x , n y ,n z ). 
This allows us to push the series further, by combining the results on various size lattices 
to cancel the contributions from particular closed loops. 

For an explicit example, consider loops of length 9. The (16,18,21) lattice has a 
minimal such loop with steps n = (3, 2, —4), the (16, 17, 21) lattice has closed loops with 
steps (1,4,-4) and (5,-1,-3), the (13,18,20) lattice has a closed loops with (2,3,-4) 
and (4, —4, 1), and finally the (14, 17, 19) system has the loops (3, 2, —4) and (5, —3, — 1). 
If we combine the coefficients as obtained from these lattices with weights (2, 1, —1, — 1) 
respectively, then all errors from the loops of length 9 cancel out. This gives the series to 
the same order as a lattice with the smallest loop having length 10, which otherwise would 
require at least 24 sites. 

This procedure extends to cancel further loops. It is straightforward to write a pro- 
gram to enumerate the closed loops on various lattices, and then solve the linear equations 
to cancel the errors from such loops. In Table I we present a list of 26 lattices and the 
relative weights for combining them to cancel all loops of length less than 14. Note that in 
this way we have reduced what would naively require a 14 3 lattice to a set of calculations 
involving a cross section of at most 24 spins. 

After cancelling the single loops as above, a potential problem arises from more com- 
plicated diagrams which wrap around the lattice simultaneously in two or more ways. This 
would correspond to flipping a set of spins which connects three of the identified sites in 
Fig. (1). In selecting our lattices for Table I we did not consider any system which had a 
loop contributing to any order for which we were extracting the series coefficient. 

It is easy to calculate the order at which these more complex loops contribute. In our 
lattice finding program we first find the three closest identified sites which do not lie on a 
single straight line. (Double loops connecting points in a line are automatically cancelled 
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at the same time as the simple loops.) Denote the minimum distances between these three 
images as dl, d2, and d3. In most cases, the minimal way of flipping a set of spins to wrap 
around these three loops produces an energy of (d — 1) x (dl + d2 + d3 — 2) — 2, where d 
is the dimension of the system. We rejected using any lattices for which this number is at 
or below the order to which we were extracting the series. 

In rare special cases this formula needs a correction. The energy can be lower if one of 
the fundamental loops has no steps in one direction. Then the two loop diagram can run 
into its periodic image, reducing the relevant order. For example, with the h = (11, 15, 18) 
lattice, the fundamental loops have n = (0, 6, —5), (3, —1, —1) and (3, 5, —6). The minimal 
energy for a set of flipped spins which connects these three images is 52 bonds, rather 
than the predicted 54 from the above formula. Needless to say, this lattice caused us 
considerable consternation. 

The utility of these cancellations depends strongly on dimension. For two dimensions 
with at most h y sites on the top row, the best solution is always a single lattice with 
h x = h y — 1. In this case the shortest extraneous loop has n = (h y , —h y + 1) with length 
2h y — 1. Note that as h y becomes large, the transfer matrix effectively grows the lattice 
along a diagonal. 

For higher dimensions, on the other hand, there are a rapidly growing number of 
interesting lattices to cancel loops between, and this method becomes particularly powerful. 
Table II includes a list of 15 lattices which give the four dimensional series through order 
50 excited bonds. Although the largest lattice here has 28 sites in the top row, the tricks 
of the next section are also more effective in four dimensions, so this is not a particularly 
difficult case. 

Note that although we have been discussing these lattices in the context of the Ising 
model, the results are more general. In particular, the combinations in Table I are valid 
for any nearest neighbor model on a simple cubic lattice. 

These methods can also be applied to other than simple cubic lattices. For example, 
to treat a body centered cubic lattice, each site has eight neighbors, so we need four 
components for h. We can merely use a four dimensional lattice-finding program modified 
to require the real closed loops of length 3 be present and not be cancelled. 
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Table I. A combination of 26 lattices which gives the three dimensional low temperature 
expansion coefficients through 54 excited bonds. The first column represents the coefficient 
with which the lattice is to be weighted and the second gives the vector h which defines 
the lattice. 



coefficient 




i 


Cl 7 23 241 


9 


C19 22 241 


1 

- 1 


Cl Q 91 941 


Q 
O 


(-\ n on o^i 


1 

_L 


Cl 8 90 931 

1 J.O ; ZiVJ ; ZjO J 


-4 


Cl 1 1 ^ 931 


_1 


CIS 21 221 


3 


Cl 6 21 221 


_2 

Zj 


Cl 8 1 9 221 




C15 19 221 




C14 19 221 


_3 


Clfi 17 221 


_3 

r > 


C5 18 211 


Q 


C8 1 7 91 1 


o 


C7 1 9 9f)l 


-6 


Cl 17 201 


6 


(16,18,19) 


-2 


(16,17,19) 


6 


(12,17,19) 


7 


(8,17,18) 


-3 


(7,16,18) 


-9 


(11,14,18) 


-7 


(8,13,18) 


2 


(9,16,17) 


3 


(1,13,15) 


4 


(12,13,14) 



Miscellaneous tricks 

During the recursive buildup of the lattice, each new count is the sum of just two 
terms, representing the two possibilities for the covered spin. Thus the arithmetic involved 
is rather trivial. On the other hand, we must store counts for all energies up to the 
maximum order desired as well as for all relevant values of the top h z spins of our helical 
lattice. In addition, the intermediate counts can become rather large numbers. Thus, the 
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Table II. A combination of 15 lattices which gives the four dimensional low temperature 
expansion coefficients through 50 excited bonds. The first column represents the coefficient 
with which the lattice is to be weighted and the second gives the vector h which defines 



the lattice. 




coefficient 


/ 7 7 7 \ 

(h x , h v , h z ) 


3 


(15,24,25,28) 


-27 


(15,21,25,28) 


14 


(13,20,25,28) 


27 


(15,20,26,27) 


16 


(11,20,26,27) 


18 


(19,20,25,27) 


2 


(11,15,25,27) 


-4 


(16,17,23,27) 


-13 


(14,17,19,27) 


-4 


(11,20,25,26) 


-6 


(15,18,25,26) 


-16 


(15,21,23,26) 


-16 


(7,20,24,25) 


28 


(14,15,23,25) 


-21 


(17,18,22,25) 



primary computational problem is storage. To substantially reduce these demands, we 
calculated the series coefficients several times, each time modulo a small different integer. 
Depending on the integers chosen, this enabled us at intermediate stages to store the 
counts in either one byte or one short integer each. As all operations are simple additions 
or multiplications, this procedure correctly gives the final coefficients modulo the given 
integers. After multiple passes using mutually prime values for these modulos, we use the 
Chinese remainder theorem to reconstruct the final series. This theorem states that if you 
know a number modulo a set of relatively prime integers, then the number is uniquely 
determined up to the product of those integers. 

As we are repeating the series calculations for several different modulos and for several 
different lattices and only combining the results at the end, this problem is particularly 
suitable for trivial parallelization. Indeed, except for the most memory intensive cases, we 
have experimented quite successfully with sending different lattice-modulo combinations to 
a farm of workstations. For this we have been using the Condor distributed batch system 
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[9]. 

Note that as we add spins, the energy of the system can only increase. This means 
that we never need counts involving more excited bonds than the order to which we are 
evaluating the series. Furthermore, while the recursive procedure is predicated on keeping 
all top rows for the lattice, this is not actually necessary if we only want the series to some 
given order. In particular, we need not store any counts for top rows which already contain 
more excited bonds than the order we are working to. To handle this, we use a simple 
subroutine that, given a possible top row, finds the next top row in numerical sequence 
with a number of excited bonds less than or equal to the working order. 

As an explicit example with the Ising model, consider the h = (17, 23, 24) lattice 
and allowing only up to 54 excited bonds. In this case we need keep only 2,778,176 of 
the possible 2 24 = 16, 777, 216 possible top rows. In four dimensions, because there are 
additional bonds which can be excited, the corresponding reductions are even larger. 

In addition to not keeping all top rows, we need not store counts with less energy 
than the minimum possible for a given top row. That is, while for the top row with all 
spins up we need to keep counts for all possible excitation energies up to the order under 
consideration, if the top row has a single flipped spin we need only keep counts of at least 
3 excited bond pairs, and so on. Finally, for the Ising case on a simple cubic lattice with 
our boundary conditions there can only be an even number of excited bonds. In this way, 
the above (17,23,24) lattice requires keeping track of 11,259,428 individual p (E,I) : or 
less than one count per possible top row. 

During the recursion, each new count is the sum of one or two of the old ones, corre- 
sponding to whether the covered spin is flipped or not, and whether for a flipped spin we 
do not already have more energy than being considered for the count in question. A simple 
way to implement this is to have two index arrays, with the elements of each representing 
the location of the old counts to be used. Having an entry in the index array out of bounds 
provides a simple way to flag those cases where only a single term goes into the sum. Once 
the geometry is established by the construction of these arrays, the program simply loops 
over the counts, making the new values the sum of two old ones pointed to by these indices. 
In this way all the complications of setting up the geometry need only be done once per 
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lattice. 

One can save additional memory by not storing the full indices, but using the fact that 
if one orders the counts first by top row numerically, and then by bonds, the respective 
indices always change by relatively small numbers in going from one count to the next. 
Thus we need only store the changes rather than the indices themselves. In the main loop 
the new indices are obtained by a simple addition to the previous ones. The index changes 
for our studies could all be stored in a single byte. 

A final trick that we have so far only used minimally is to invert a partially grown 
lattice on itself. The idea is that given the counts for all possible top layers, we can then 
obtain the counts for a lattice roughtly twice as long with all possible specified layers in 
the middle. Calling this count Pd(E, I), we have 

p d (E, I) = P( E ^ T )P( E ^ ^(E, E x + E 2 - d{I)) (11) 

Ei,E2 

where I r has the bits of I in reversed order (because the lattice has been flipped upside 
down) and d(I) represents the excited bonds inside the middle layer. The latter is removed 
to prevent double counting. This technique can provide information on correlation func- 
tions in this middle layer. As all states are known explicitly, any such correlation function 
can be obtained exactly with no significant additional drains on computer time or memory. 
As a simple example, this provides an alternative method for obtaining the magnetization 
series to that discussed in the next section. 

Other observables 

So far we have been discussing the direct low temperature expansion for the partition 
function or, equivalently, the average energy or the specific heat. The method easily 
extends to other observables by generalizing the counts. For example, consider applying a 
magnetic field by generalizing the partition function to 

z = £y/»*-ffI>« (12) 

Derivatives with respect to the applied field give us a procedure to compute the 
magnetization 

M = (a i ) = -±-^\ogZ (13) 
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and the magnetic susceptibility 

x = Wh m (14 > 

For general H one can expand observables simultaneously in u and A = e~ 2H . In 
Ref. [6] this possibility was discussed in terms of generalizing the counts P(E) to the 
two indexed count P(E 7 S) : representing the number of states of a given bond energy 
and number of flipped spins S. The recursion relations for these counts are completely 
analagous to those for P(E). The double series for the magnetization was presented up to 
order 42 excited bonds in Ref. [6]. 

One difficulty with this approach is the increased memory required for storing counts 
for all magnetizations as well as energies. If one is only inte rested in the magnetic proper- 
ties in the zero field limit, one can store considerably less. In particular, consider moments 
of the magnetization, from which quantities such as the susceptibility are easily extracted. 
It is convenient to define new quantities 

P k (E) = J2S k P(E,S). (15) 

s 

With this definition, Pq(E) is simply the original count P{E). The zero field magnetization 
is easily found from 



1 2 E E Pi(E)e- pE 



M = l-2 ^ • (16) 

Finally, from P 2 we can obtain the magnetic susceptibility 

x = 4 ( E B (^)-fi(S)')e-»« ). 

The advantage of working with these moments is that they themselves satisfy simple 
recursion relations. To derive them, consider the generalization of Eq. (4) 

p'(E, S, I) = $>(£- I%S- A.(I'), I')- (18) 
r 

Here p(E, S, I) is the number of states of energy E, with S flipped spins, and with lattice 
top row specified by /, and p' is the same quantity on the new lattice obtained after adding 
the new spin. We denote by A S (I') the change in the number of flipped spins; that is, 
A s = 1 if the new spin is flipped (the relevant bit of V = 1) and A s = otherwise. 
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Now define the moments, 



Pk (E,I) = Y,S k P (E,S,I). 



(19) 



s 



Taking moments of Eq. (18) now gives the recursion relations for the p k 



Po(EJ) 



J2po(E-A,I') 



(20) 



i' 



£>i(£ - A, I') + A s (I)p (E - A, I') 



(21) 



p' 2 (E, I) = Y J P2{E- A, I') + 2A s (I) Pl (E - A, I') + A 2 s (I)p (E - A, /')• (22) 



The first of these relations is just our original recursion, and the others enable us to 
calculate the magnetization and susceptibility with the addition of only two new counts. 

It is straightforward to derive the analogous counting schemes for n-point susceptibil- 
ities and their various spatial moments, like the second moment of 2-point susceptibility 
H2 =< x 2 a x ao >. In later case however there are some conceptual difficulties connected 
to the ambiguity of the definition of the coordinate on the helical lattice. Some more work 
needs to be devoted to this problem. 

Strong coupling 

We now turn to the application of the counting methods to the strong coupling series. 
In this section we describe the procedure for the three dimensional Ising model, although 
again it is easily generalized. As before we consider spins Si on the lattice sites i and 
taking the values ±1. The partition function of Eq.(2) can be trivially rewritten 



where the product is over all lattice links and Ni is the number of links in the system. 
For the strong coupling series we consider small (3 and expand the above sum in powers of 
tanh(/3). Each term involves a set of selected bonds which each give a power of tanh(/3). 
Having selected a set of bonds, we can then perform the sum over the spins. If any site 




(23) 
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has an odd number of selected bonds eminating from it, the sum will vanish. Otherwise 
the sum over any given spin gives a factor of two. Thus we conclude 

Z = 2 N (^—) N > £ iV(fc)(tanh(/3)) fc (24) 
1 k 

Here N(k) represents the number of possible ways to select k bonds in such a manner that 
each site is the end of an even number of selected bonds. We adapt our counting methods 
to evaluate these numbers N(k). 

As before we maintain information on the top layer of our lattice while adding new 
sites one at a time. Here, however, rather than the values of the spins themselves on the 
top layer, we keep information on the selected bonds ending there. In particular, because 
we want to allow future bonds to extend above the top row, we relax the constraint that 
an even number of bonds end on the top sites. Thus, we keep a count N(k, I) where I 
now stores in its set bits those sites with an odd number of bonds coming into them from 
previous sites. We refer to sites with an odd number of incoming bonds as having "loose 
ends" or "dangling bonds." On adding a new site, we have the basic recursion relation 

JV'(M) = X>(*-A(I, !'),/') (25) 
r 

where A (J, /') represents the number of selected bonds attached to the new spin and /' is 
related to I via changes in those bits representing sites attached to the new one. 

In three dimensions, for any given (k, I) there will be four terms in the above sum 
over I'. This represents a factor of two for whether the new x bond is selected times a 
factor of two for whether the new y bond is chosen. Whether the corresponding z bond is 
chosen or not is determined by the corresponding bit of I which determines if an even or 
odd number of bonds are selected. 

An immediate factor of two in memory is saved because each bond has two ends. This 
means that if no bonds enter from outside below the lattice, the top layer must have an 
even number of loose ends. Any top layers with an odd number of loose ends need never be 
kept. In practice, instead of looping over all given any integers I representing the dangling 
bonds from an allowed configuration, we need only loop over the right h z — l bits of / and 
can determine the allowed leftmost bit by parity considerations. 
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We work with generalized helical lattices as before. For simplicity in initialization, we 
set all counts to zero except for 1 = 0, representing no dangling bonds. This may seem a 
bit peculiar because we do not allow loops to enter and travel through the bottom layer. 
It is, however, simple to implement and boundary conditions in the longitudinal direction 
are irrelevant if we grow the lattice long enough. 

On a single helical lattice, the strong coupling series will be correct to the order of the 
first chain of bonds which wraps around one of the artificial closed loops discussed earlier. 
The double loop criterion is somewhat different now; here it is only the total length of a 
loop which wraps around two directions that matters. Rejecting lattices with such double 
loops, we can perform the same cancellation between lattices as in the low temperature 
series. 

The strong coupling series can be extended significantly by using the fact that all 
valid loops of links on an infinite lattice will have an even number of selected bonds in 
any of the coordinate directions. We use this fact by calculating the counts several times, 
but including extra minus signs when adding bonds in various directions. For example, 
if we first find the series giving every x bond a weight of -1, we can then add the result 
without this extra sign and any artificial diagram involving an odd number of x bonds will 
cancel out. Thus we need not worry about any finite size effects involving an odd number 
of steps in the x direction. Repeating the procedure 8 times for all combinations of minus 
signs for the three possible directions, we can ignore any extraneous closed loops with an 
odd length along any dimension. Similarly, any double loops with an odd number of steps 
in any direction can also be ignored. Without this trick the order to which the strong 
coupling series can be found is rather uninteresting. 

With these tricks, we have found the series through 20 selected bonds from the com- 
bination of lattices given in Table (III). As the lattice size goes to infinity, we write the 
free energy in the form 

F = ^ = log 2 + ^ log(i±|^!) + £ /fc (tanh(/3)) fc (26) 

k 

where Ns and Nl denote the number of sites and links, respectively. To extract the 
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Table III. A combination of 11 lattices which gives the three dimensional strong coupling 
series through order 20. The first column represents the coefficient with which the lattice 
is to be weighted and the second gives the vector h which defines the lattice. 

coefficient (h x ,h y ,h z ) 

4 (4,15,16) 

2 (12,13,16) 
-4 (4,13,16) 
-6 (11,12,16) 
4 (7,12,16) 
6 (11,14,15) 
-4 (7,13,15) 

3 (5,11,15) 
3 (9,13,14) 
-3 (4,12,13) 
-3 (8,11,13) 



Table IV. The coefficients for the strong coupling series for the three dimensional Ising 
model through order 20. The are defined in the text. 



k kf k 



2 

4 12 

6 132 

8 1,500 

10 19,800 

12 288,528 

14 4,468,380 

16 72,236,124 

18 1,206,062,448 

20 20,649,134,532 



coefficients it is somewhat easier to work with the analog of an expectation, 

As for the low temperature series, we extract the contribution per spin by comparing the 
counts before and after adding the last spin. Since they are just combinations of integer 
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counts, the products kfk themselves are always integers, while the //- are not in general. 
We tabulated these numbers through order 20 in Table (IV) . These numbers are not new; 
for example they plus one additional term follow from the results in Ref. [3]. 

Potts models 

As we mentioned earlier, the application of the counting techniques to the low tem- 
perature series expansions is very easily generalizable to any discrete system with nearest 
neighbor interaction. To illustrate this, consider the q-state Potts model, defined by the 
interaction of the form 

£ = £[1-<WJ (28) 

(ij) 

where Oi is a site-defined field that takes q possible values. The sum is taken over all 
nearest neighbor pairs of spins with 6 being the Kronecker symbol. 
Writing the partition function in the form 

dN 

Z=J2 P(E)u E (29) 

E=0 

with d being the spatial dimension and u = e - ^, one can follow essentially the same steps 
we outlined in the discussion of the Ising model. Namely, the application of recursive 
counting using helical lattices and Chinese arithmetic comes through with no change at 
all. The differences are of a technical nature only, the conceptual ones. 

Working with h spins on a helix, the maximum number of configurations of the top 
layer is q . Since a single bit is no longer sufficient to keep the state of the individual 
spin, it would be more complicated to code the state of the top layer in a single word. 
Instead, we use several words to represent each top layer configuration. For example, for 
the q = 3 calculations we used two words per configuraion while for q = 8 three words 
were required. It is also clear that now the analogue of Eq. (4) has q terms, corresponding 
to the q different possible values of the added spin. 

We have computed the low temperature expansions for the energy, magnetization and 
susceptibility for the q = 3 model in d = 2 and d = 3 and the q = 8 model in d = 2. The 
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resulting series have been extensively discussed and analyzed in [7] and we do not repeat 
them here. 

Results and analysis 

Using these methods with the simple Ising model, we obtained the series for the 
average energy per bond given in Table (V). In Table (VI) we give the series through 
order 54 excited bonds for the magnetization and the magnetic susceptibility of the three 
dimensional model. 

Table V. The low temperature expansion coefficients for the average energy per unit volume 
for the three, four, and five dimensional Ising model on a simple cubic lattice. 



% 


p- (1 rH 


f>- (A rH 
















2 











4 











6 


12 








8 





16 





10 


60 





20 


12 


-84 








14 


420 


112 





16 


-1,056 


-144 





18 


3,756 





180 


20 


-11,220 


1,120 


-220 


22 


37,356 


-2,816 





24 


-118,164 


2,032 





26 


389,220 


11,856 


2,340 


28 


-1,261,932 


-46,704 


-5,600 


30 


4,163,592 


66,960 


3,320 


32 


-13,680,288 


94,576 


640 


34 


45,339,000 


-707,472 


32,980 


36 


-150,244,860 


1,545,120 


-122,220 


38 


500,333,916 


-148,656 


145,540 


40 


-1,668,189,060 


-9,522,864 


-31,420 


42 


5,579,763,432 


30,130,576 


454,860 


44 


-18,692,075,820 


-30,299,808 


-2,483,360 


46 


62,762,602,860 


-104,198,096 


4,560,440 


48 


-211,062,133,044 


520,429,776 


-2,922,240 


50 


711,052,107,060 


-918,744,400 


6,717,220 


52 


-2,398,859,016,684 






54 


8,104,930,537,260 
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Table VI. The low temperature expansion coefficients for the average magnetization and 
magnetic susceptibility for the three dimensional Ising model on a simple cubic lattice. 
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Mi 


Xi 
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n 





n 

U 
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A 



u 


n 

w 


6 

U 




Zj 


1 


8 
o 
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1 2 


1 2 


1 9 


14 


14 


14 


-Q0 


1 35 

-LOU 


1 fi 


1 Q2 

J. C/Zj 


Zj / U 




-7Q2 

1 t7 Zj 


1 520 

-L • (J Zj \J 


90 

Zj W 


2 148 

Zj . X t:0 


-4 05fi 


99 


-7 71 


1 7 778 


24 

Zj^ 


2S 262 


-54 SQ2 


26 

Zj \) 


-7Q 51 2 


21 3 522 

Zj _L kj i (J Zj Zj 


28 

Zj O 


252 054 

Zj tj Zj . U(Jt: 


-700 362 


30 


-84fi fi28 


2 fiOl fi74 

Zj . Uu J. ; U f L ± 


32 

OZv 


2 753 520 


-8 83fi 81 2 


34 


-9 205 800 


31 925 04fi 


36 


SO 371 124 


-110 323 056 


38 


-101 585 544 


393 008 71 2 

Oi/0 ; VJWO ; / -LZj 


40 


338 095 596 


-1 369 533 048 


42 


-1,133,491,188 


4,844,047,090 


44 


3,794,908,752 


-16,947,396,000 


46 


-12,758,932,158 


59,723,296,431 


48 


42,903,505,030 


-209,328,634,116 


50 


-144,655,483,440 


736,260,986,208 


52 


488,092,130,664 


-2,582,605,180,212 


54 


-1,650,000,819,068 


9,074,182,912,884 



We will now give a brief analysis of our series to get results on the critical temperature 
and exponents. In the usual Dlog Pade (DIP) analysis [10], given a series expansion for 
F{u) to iV-th order, F^{u) = 1 + J^Li f% u% i ( we wu l use the simplification that one can 
always normalize the series so that the constant term is unity), one computes coefficients 
for polynomials Ql(u) = Yl^=Q ( li u% and Rm{u) = 1 + Yl^iLi 7 "^ '■> which satisfy, 

Ql(u)/R m (u) = F' N {u)/F N {u) (30) 

to 0(u N ) with L + M = N - 1. 
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The position u c of a singularity of F of the form F ~ A/\u — u c \^ will be approximated 
by the zeros of Rm- The exponent ( is estimated by £ = —Q l(u c ) / R'm( u c) ■ 

In addition to DIP, we will also use the method of inhomogeneous differential approx- 
imants (IDA) introduced by Fisher and Au-Yang [11] (see also [12]). These are useful in 
handling singularities of the form, 

F(u) = , Ai - U \ t +B(u) (31) 
\u — u c \^ 

where A and B are analytic in u. In this method, one computes coefficients for polynomials 
Ql(u), Rm(u) and Sj(u) which satisfy, 

FnQl + Sj = F N R M (32) 

to order N, with L + M + J = N — 2. Here the subscripts represent the highest order 
present in the given polynomials. 

Note that for Sj = one gets the usual Dlog Pade ratio from Ql/Rm- It is easy to 
see that potential critical points u c are the zeros of Rm and for each of these, the exponent 
C is again estimated as ( = —Ql(u c )/R'm( u c)- 

In Fig. 2 we show the critical point u c and the exponent (3 obtained from the magne- 
tization series in Table (VI) using DIP and IDA analysis. The different points represent 
different values for the orders of the various polynomials used in the analysis. The crit- 
ical value u c is most accurately known from Monte Carlo calculations and is shown as a 
vertical line. The thickness of this line represents the accuracy with which this number is 
known. As is clear, the DIP gives a u c quite far from the Monte Carlo value and the IDA 
results tend to scatter around it. If however, we take u c as given, we can calculate (3 quite 
accurately from our results by interpolating the IDA results to the known u c . Fitting a 
straight line to the points near u c in Fig. (2) gives (3 = 0.289(1) where the errors are only 
the errors on the fit and ignore possible systematic effects from the unknown higher order 
terms in the series. 

Similarly, the susceptibility series gives the results in Fig. 3. Now the results for the 
DIP and the IDA are mixed in together, in contrast to the magnetization series results 
(see above). Once again, the critical point is more accurately obtained from Monte Carlo 
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data, and, assuming that value, we find the exponent 7 = 1.281(1) from straight line fits 
to the data in the critical region. 

Finally, we turn to the heat capacity, C v , series. The results are shown in Fig. (4). 
The DIP lie deceptively on a straight line which, if extrapolated to the Monte Carlo value 
for u c gives the wrong answer near a = 0.2 as was already noticed in our earlier paper [13]. 
The IDA results on the other hand, are steeply varying in the region of u c . If we take the 
three points from the J = IDA data and fit them to a line, we get a = 0.128(1). 

Concluding remarks 

The method presented here should easily generalize to other discrete systems. The 
helical lattices used, as well as the combinations to cancel out finite size errors, are in- 
dependent of the specific model. It is straightforward to introduce additional couplings, 
although this will increase memory needs. Some interesting possibilities for futher explo- 
ration are gauge and coupled gauge-spin models in various dimensions. Changing boundary 
conditions should enable the study of interface properties. In Ref. [14] similar recursive 
methods were suggested as a means to study many fermion systems. A particularly chal- 
lenging problem is the extension of these ideas to theories with continuous spins. Some 
work along these lines for gauge theories appears in Ref. [15] 
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Figure captions 

Fig. 1. Visualizing the helical lattice. All lattice points lying in the plane orthogonal 
to the vector h are to be identified. 

Fig. 2. The critical coupling and exponent f3 obtained from the magnetization series 
for the three dimensional Ising model. 

Fig. 3. The critical coupling and exponent 7 obtained from the magnetic susceptibility 
series for the three dimensional Ising model. 

Fig. 4. The critical coupling and exponent a obtained from the energy series for the 
three dimensional Ising model. 
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